課題の回答例

【課題】限界集落の可視化

e-Stat(統計GIS)の平成22年国勢調査(少地域)から適当な地域を選び(出身地でも国立市でも、鯖江市以外なら良い)、「地域別高齢者(65歳以上)人口比率」がみえるように可視化してください。

可視化例:

if(!require(maptools)){
  install.packages("maptools")
  require(maptools)
}
if(!require(ggmap)){
  install.packages("ggmap")
  require(ggmap)
}
# ggplot2はggmap 読み込み時に自動的に読み込まれるが念のため
if(!require(ggplot2)){
  install.packages("ggplot2")
  require(ggplot2)
}

#日本語フォント設定(Win/Mac両対応)
if(.Platform$OS.type=="windows")
  windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
  quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))

#福井県鯖江市、平成22年国勢調査(少地域)2010/10/1
# shapefileの読み込み
shp<-maptools::readShapePoly("A002005212010DDSWC18207/h22ka18207.shp")
# data.frameへの変換
dat <- broom::tidy(shp)
# 行番号を新規列項目として追加
shp@data$id <- rownames(shp@data)
# 「dat」と「shp@data」とを列「id」をキーに結合
shp.df <- dplyr::inner_join(dat, shp@data,by = "id")
# 列「KEY_CODE」を数値に変換
shp.df$KEY_CODE<-as.numeric(as.character(shp.df$KEY_CODE))
# 統計データの読み込み
csv<-read.csv("tblT000573C18207.txt",stringsAsFactors = FALSE,fileEncoding = "SJIS")
# 内容確認
head(csv)
##     KEY_CODE HYOSYO CITYNAME       NAME HTKSYORI HTKSAKI GASSAN
## 1         NA     NA                                   NA     NA
## 2 1.8207e+04      1   鯖江市                          NA     NA
## 3 1.8207e+08      2   鯖江市       本町               NA     NA
## 4 1.8207e+10      3   鯖江市 本町1丁目               NA     NA
## 5 1.8207e+10      3   鯖江市 本町2丁目               NA     NA
## 6 1.8207e+10      3   鯖江市 本町3丁目               NA     NA
##               T000573001   T000573002   T000573003       T000573004
## 1 総数、年齢「不詳」含む 総数0〜4歳 総数5〜9歳 総数10〜14歳
## 2                  67450         3344         3604             3595
## 3                    871           26           35               36
## 4                    184           10            6               10
## 5                    218            2            7               13
## 6                    121            5            4                6
##         T000573005       T000573006       T000573007       T000573008
## 1 総数15〜19歳 総数20〜24歳 総数25〜29歳 総数30〜34歳
## 2             3335             2768             3603             4487
## 3               30               20               26               50
## 4                4                5               10               11
## 5               10                7                7                7
## 6                2                4                3               16
##         T000573009       T000573010       T000573011       T000573012
## 1 総数35〜39歳 総数40〜44歳 総数45〜49歳 総数50〜54歳
## 2             5303             4291             4038             3811
## 3               45               30               36               57
## 4               10                7                4               10
## 5               13                7               11               12
## 6                6                3                7               11
##         T000573013       T000573014       T000573015       T000573016
## 1 総数55〜59歳 総数60〜64歳 総数65〜69歳 総数70〜74歳
## 2             4421             5373             4126             3470
## 3               52               81               65               63
## 4               14               20               16               13
## 5               11               12               17               17
## 6                8               14                8                7
##       T000573017       T000573018     T000573019     T000573020
## 1 総数15歳未満 総数15〜64歳 総数65歳以上 総数75歳以上
## 2          10543            41430          15371           7775
## 3             97              427            347            219
## 4             26               95             63             34
## 5             22               97             99             65
## 6             15               74             32             17
##                   T000573021 T000573022 T000573023     T000573024
## 1 男の総数、年齢「不詳」含む 男0〜4歳 男5〜9歳 男10〜14歳
## 2                      32507       1686       1830           1846
## 3                        386         14         21             14
## 4                         84          4          2              4
## 5                         93          1          5              5
## 6                         53          3          1              1
##       T000573025     T000573026     T000573027     T000573028
## 1 男15〜19歳 男20〜24歳 男25〜29歳 男30〜34歳
## 2           1768           1292           1761           2211
## 3             15              9             14             25
## 4              3              2              6              6
## 5              5              2              4              3
## 6              -              3              1              9
##       T000573029     T000573030     T000573031     T000573032
## 1 男35〜39歳 男40〜44歳 男45〜49歳 男50〜54歳
## 2           2723           2141           1992           1867
## 3             22             14             17             24
## 4              7              1              4              1
## 5              6              4              4              6
## 6              2              2              3              6
##       T000573033     T000573034     T000573035     T000573036   T000573037
## 1 男55〜59歳 男60〜64歳 男65〜69歳 男70〜74歳 男15歳未満
## 2           2124           2595           2090           1598         5362
## 3             24             37             35             27           49
## 4              6             11              9              6           10
## 5              6              4             10              7           11
## 6              2              6              3              4            5
##       T000573038   T000573039   T000573040                 T000573041
## 1 男15〜64歳 男65歳以上 男75歳以上 女の総数、年齢「不詳」含む
## 2          20474         6597         2909                      34943
## 3            201          136           74                        485
## 4             47           27           12                        100
## 5             44           38           21                        125
## 6             34           14            7                         68
##   T000573042 T000573043     T000573044     T000573045     T000573046
## 1 女0〜4歳 女5〜9歳 女10〜14歳 女15〜19歳 女20〜24歳
## 2       1658       1774           1749           1567           1476
## 3         12         14             22             15             11
## 4          6          4              6              1              3
## 5          1          2              8              5              5
## 6          2          3              5              2              1
##       T000573047     T000573048     T000573049     T000573050
## 1 女25〜29歳 女30〜34歳 女35〜39歳 女40〜44歳
## 2           1842           2276           2580           2150
## 3             12             25             23             16
## 4              4              5              3              6
## 5              3              4              7              3
## 6              2              7              4              1
##       T000573051     T000573052     T000573053     T000573054
## 1 女45〜49歳 女50〜54歳 女55〜59歳 女60〜64歳
## 2           2046           1944           2297           2778
## 3             19             33             28             44
## 4              -              9              8              9
## 5              7              6              5              8
## 6              4              5              6              8
##       T000573055     T000573056   T000573057     T000573058   T000573059
## 1 女65〜69歳 女70〜74歳 女15歳未満 女15〜64歳 女65歳以上
## 2           2036           1872         5181          20956         8774
## 3             30             36           48            226          211
## 4              7              7           16             48           36
## 5              7             10           11             53           61
## 6              5              3           10             40           18
##     T000573060
## 1 女75歳以上
## 2         4866
## 3          145
## 4           22
## 5           44
## 6           10
# 先頭行の削除
csv <- csv[-1,]
dat <- data.frame(KEY_CODE=csv$KEY_CODE,TOTAL=csv$T000573019)
shp.df2 <- dplyr::inner_join(shp.df,dat,by="KEY_CODE")
shp.df2$TOTAL <- as.numeric(as.character(shp.df2$TOTAL))
shp.df2$RATE <- shp.df2$TOTAL/shp.df2$JINKO
# Inf(無限)となっていたものを0に置換(0除算結果の修正)
shp.df2$RATE[is.infinite(shp.df2$RATE)] <- 0
#NA(欠損値)を0に置換
shp.df2$RATE[is.na(shp.df2$RATE)] <- 0

#ラベル表示の準備
#文字コード変換(SJIS > UTF8).Windows環境の場合は不要.
shp.df2$MOJI <- iconv(as.character(shp.df2$MOJI),from="SJIS",to="UTF8")
#ラベル座標とラベルの設定
shp.text <- data.frame(shp.df2$X_CODE,shp.df2$Y_CODE,shp.df2$MOJI)
#重複を除去する
shp.text <- unique(shp.text)
names(shp.text) <- c("lon","lat","name")
shp.text$name <- as.character(shp.text$name)
#GoogleMapを取得する
map <- ggmap::get_map(location=c(136.2302503,35.9633669),zoom=12,maptype="roadmap")
#GoogleMapを描画する(ggplot2の初期化も実施している)
gm <- ggmap::ggmap(map)
#市町村の領域を描画.colorで境界線の色を指定.alphaで透明度を指定
g <- gm + ggplot2::geom_polygon(data = shp.df2, aes(x = long, y = lat, group = id,fill=RATE),color="gray80",alpha=0.7)
#塗り色を変更する
g <- g+ ggplot2::scale_fill_gradientn(colours=rev(heat.colors(100)))
#タイトルを指定する.
g <- g +ggplot2::ggtitle("福井県鯖江市の高齢者比率")+ggplot2::theme_gray(base_family = "yugo")
plot(g)

【発展課題】市営駐車場の可視化

以下の3つのデータを組合せた可視化にトライして下さい.

完成イメージは以下のとおりです。

if(!require(maptools)){
  install.packages("maptools")
  require(maptools)
}
if(!require(ggmap)){
  install.packages("ggmap")
  require(ggmap)
}
# ggplot2はggmap 読み込み時に自動的に読み込まれるが念のため
if(!require(ggplot2)){
  install.packages("ggplot2")
  require(ggplot2)
}
#XML処理に必要
if(!require(XML)){
  install.packages("XML")
  require(XML)
}
#shapefileを読み込む
shp<-maptools::readShapePoly("A002005212010DDSWC18207/h22ka18207.shp")
#data.frameに変換
dat <- broom::tidy(shp)
# 行番号を新規列項目として追加
shp@data$id <- rownames(shp@data)
# 「dat」と「shp@data」とを列「id」をキーに結合
shp.df <- dplyr::inner_join(dat, shp@data,by = "id")
# 列「KEY_CODE」を数値に変換
shp.df$KEY_CODE<-as.numeric(as.character(shp.df$KEY_CODE))
# 統計データの読み込み
csv<-read.csv("tblT000573C18207.txt",stringsAsFactors = FALSE,fileEncoding = "SJIS")
csv <- csv[-1,]
# 「csv$KEY_CODE,csv$T000573019」からなるdata.frameを作る
dat <- data.frame(KEY_CODE=csv$KEY_CODE,TOTAL=csv$T000573019)
# 「dat」と「shp.df」とを列「KEY_CODE」をキーに結合
shp.df2 <- dplyr::inner_join(shp.df,dat,by="KEY_CODE")
# 数値に変換(計算のため)
shp.df2$TOTAL <- as.numeric(shp.df2$TOTAL)
# 割合の計算
shp.df2$RATE <- shp.df2$TOTAL/shp.df2$JINKO
# エラー処理
shp.df2$RATE[is.infinite(shp.df2$RATE)] <- 0
shp.df2$RATE[is.na(shp.df2$RATE)] <- 0
#文字コード変換(SJIS > UTF8)
shp.df2$MOJI <- iconv(as.character(shp.df2$MOJI),from="SJIS",to="UTF8")
shp.text <- data.frame(shp.df2$X_CODE,shp.df2$Y_CODE,shp.df2$MOJI)
# 重複除去
shp.text <- unique(shp.text)
names(shp.text) <- c("lon","lat","name")
shp.text$name <- as.character(shp.text$name)

if(!file.exists("parking.xml")){
  download.file(url="http://www3.city.sabae.fukui.jp/xml/parking/parking.xml",destfile="parking.xml")
} 
#XMLの読み込み
doc<-XML::xmlParse("parking.xml",encoding="UTF8",useInternalNodes=TRUE)
#要素「parkinng1」を数える
len <- length(XML::getNodeSet(doc,"//parkinng1"))
#空のdata.frameを作る
dat <- data.frame()
#要素「parkinng1」の分繰り返す
for(i in 1:len){
  # 要素値を取得する
  label <- XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/name"))[[1]])
  latitude <- as.numeric(XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/latitude"))[[1]]))
  longitude <- as.numeric(XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/longlatitude"))[[1]]))
  size <- as.numeric(XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/max"))[[1]]))
  # 取得した要素値からdata.frameを作る
  item <- data.frame(label,latitude,longitude,size)
  dat <- rbind(dat,item)
}
# Googlemapを取得する
map <- ggmap::get_map(location=c(136.2302503,35.9633669),zoom=12,maptype="roadmap")
g <- ggmap::ggmap(map)
g <- g + ggplot2::geom_polygon(data = shp.df2, aes(x = long, y = lat, group = id,fill=JINKO),color="gray80",alpha=0.7)
g <- g + ggplot2::geom_point(data=dat,aes(x=longitude,y=latitude,size=size),color="blue")
g <- g+ ggplot2::scale_fill_gradientn(colours=rev(heat.colors(100)))
plot(g)

はじめに

オープンデータについてもう一度振り返ります. オープンデータとは「オープンなデータ」のことです. Open Definitionでは「誰もが自由に利用、再利用、再配布可能なデータ」と定義しています.

これまで演習では「オープンデータの評価指標「5 Star OPEN DATA(以下図参照)」に沿ったオープンデータの可視化に挑戦してきました.とくに三ツ星までの「オープンライセンス(Open License;OL)」「マシンリーダブル(Machine Readable;RE)」「オープンフォーマット(Open Format;OF)」の3要素が重要になります.

まず,そのデータがオープンであることを示すには,Creative Commons([https://creativecommons.org/](https://creativecommons.org/) などのオープンライセンスを用いてオープンであることを明示する必要があります.

つぎに,データとして公開するのであれば, マシンリーダブル,すなわちコンピュータで処理しやすい形で提供する必要があります. PDFで公開されたデータより,Excelで公開されたデータの方が便利ですよね.

そして,データフォーマットもプロプライエタリ(データ形式等の仕様が非公開)なフォーマットよりは, オープンフォーマットの方が,処理可能性が広がり,誰もが利用できるデータとなるのです.

また,どんなデータでもオープンデータとなり得るのですが, その中でも価値あるデータと位置付けられているものがあります. G8オープンデータ憲章には「ハイバリューデータ・セット」として 利用価値の高いデータ種類がまとめられています.

データカテゴリ データセット例
企業 企業/事業者の登記情報
犯罪と司法 犯罪統計、安全
地球観測 気象/天候、農業、林業、水産業及び狩猟
教育 学校のリスト; 学校の業績、デジタルスキル
エネルギーと環境 汚染レベル、エネルギー消費
財政と契約 取引支出、賃貸契約、入札需要、将来入札、地方公共団体予算、国家予算(執行計画及び決算)
地理空間 地勢、郵便番号、国の地図、地方(機関)の地図
国際開発 援助、食料安全保障、資源採掘、土地
政府の説明責任と民主主義 政府の窓口,選挙結果、法律及び規則、給与(給与水準)、福利厚生/贈与
健康 処方箋データ、実績データ
科学と研究 ゲノムデータ、研究及び教育活動、実験結果
統計 全国統計、国勢調査、社会基盤、財産、スキル
社会的流動性と福祉 居住、健康保険及び失業手当給付金
交通とインフラ 公共交通機関の時刻表、アクセスポイント、 ブロードバンドの普及度

ref. http://www.kantei.go.jp/jp/singi/it2/densi/dai4/sankou8.pdf

これまでこうした類のデータは専門家だけのものだったのですが, オープンデータやオープンサイエンスといったアイデアの浸透により,だれもが利用できるようになってきました.演習では「地球観測」「地理空間」「科学と研究」「統計」といったデータの可視化に挑戦してきましたが,なんとなくオープンデータの可能性については感じとってもらえたのではないでしょうか.

本日はRを使ったオープンデータの可視化の最終回となります. これまで以上に高度な内容が含まれていますが,挑戦する気概で演習に取り組んで下さい. 様々なデータの可視化については課題のヒントとしてご利用下さい.

Google Earthを使った可視化

これまでは ggplot2 や ggmap を使った可視化を試してきましたが, plotKML((http://plotkml.r-forge.r-project.org/)[http://plotkml.r-forge.r-project.org/])を用いると、Google Earthを用いた可視化が手軽にできます.

plotKML はその名のとおり,Google Earthで採用されているKML(Keyhole Markup Language)としてプロット結果を出力するライブラリです. 本演習では簡単な使い方を紹介します.興味がある方はチュートリアル(英語)リファレンス(英語)等をご参照ください.

shapefileを使った可視化

さっそくですが,plotKMLの実行例です.統計GISからダウンロードしたshapefileを使って地域別人口を可視化します.

if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
# KMLファイル用のライブラリ
if(!require(plotKML)){
  install.packages("plotKML")
  library(plotKML)
}
# shapefileを読み込む
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
# 投影・座標系情報を設定する。
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
# JINKOの値を高さに用いる
z<-shp@data$JINKO
z[z==0]<-1
# KML形式でplot
plotKML(shp,altitude=z)

コードは単純ですが実行してみると, Google Earthが起動し,プロット結果を表示します.

加えて地名を表示してみます。

if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
# KMLファイル用のライブラリ
if(!require(plotKML)){
  install.packages("plotKML")
  library(plotKML)
}
# shapefileを読み込む
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
# 投影・座標系情報を設定する。
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
# JINKOの値を高さに用いる
z<-shp@data$JINKO
z[z==0]<-1
#ラベル表示の準備
#文字コード変換(SJIS > UTF8)
shp@data$MOJI <- iconv(as.character(shp@data$MOJI),from="SJIS",to="UTF8")
# KML形式でplot
plotKML(shp["MOJI"],altitude=z)

散布図っぽく人口を可視化

今度はGoogle Earthで散布図もどきを描いてみましょう.

先ほどと同じように統計GISのshapefileを使います.

if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
if(!require(plotKML)){
  install.packages("plotKML")
  library(plotKML)
}
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
#必要な情報だけをもったdata.frameを作成する。
data <- data.frame(x=shp@data$X_CODE,y=shp@data$Y_CODE,id=shp@data$KEY_CODE,label=iconv(shp@data$MOJI,to="UTF8",from="SJIS"),jinko=shp@data$JINKO)
#座標を設定する。
coordinates(data) <- ~x+y
# 投影・座標系情報を設定する。
proj4string(data)<-CRS("+proj=longlat +datum=WGS84")
#KMLファイルを出力する。点(アイコン)のサイズは人口。
kml(obj=data,folder.name="",file.name="data.kml",shape="http://maps.google.com/mapfiles/kml/pal2/icon18.png",colour="#ff0000",alpha=1.0,size=jinko,labels=label)

このように緯度経度情報を含んだdata.frameさえ用意すれば, ggplot2 と同じようにプロットすることができます.

shapeファイルに使えるアイコンセットは以下のサイトで確認することができます.

http://kml4earth.appspot.com/icons.html

data.frameを用意せず,kml関数に直接引数を渡すといった使い方も可能です.

if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
if(!require(plotKML)){
  install.packages("plotKML")
  library(plotKML)
}
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
z<-shp@data$JINKO
z[z==0]<-1
kml(obj=shp,folder.name="",file.name="data2.kml",labels=iconv(shp@data$MOJI,to="UTF8",from="SJIS"),altitude=z,colour="#ff00ff",alpha=0.75,plot.lab=TRUE)

高齢者人口比率の可視化

最後に高齢者人口比率をplotKMLで可視化します.

if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
if(!require(plotKML)){
  install.packages("plotKML")
  library(plotKML)
}
if(!require(plyr)){
  install.packages("plyr")
  library(plyr)
}
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212")
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
csv<-read.csv("tblT000573C13212.txt",stringsAsFactors = FALSE,fileEncoding = "SJIS")
csv <- csv[-1,]
dat <- data.frame(KEY_CODE=csv$KEY_CODE,TOTAL=csv$T000573019)
shp@data <- join(shp@data,dat,by="KEY_CODE")
shp@data$TOTAL <- as.numeric(as.character(shp@data$TOTAL))
shp@data$RATE <- shp@data$TOTAL/shp@data$JINKO
# Inf(無限)となっていたものを0に置換(0除算結果の修正)
shp@data$RATE[is.infinite(shp@data$RATE)] <- 0
#NA(欠損値)を0に置換
shp@data$RATE[is.na(shp@data$RATE)] <- 0
#ラベル表示の準備
#文字コード変換(SJIS > UTF8)
shp@data$MOJI <- iconv(as.character(shp@data$MOJI),from="SJIS",to="UTF8")
#カラーパレットの準備
pal<-rev(heat.colors(100))
#確認
print(pal)
#カラーパレットを使って色を指定
kml(obj=shp,folder.name="",file.name="data3.kml",colour_scale=pal,colour=shp@data$RATE,labels=shp@data$MOJI,altitude=100,alpha=0.78,plot.lab=TRUE)

ggplot2のときと可視化内容は同じですが、 Google Earthなのでインタラクティブに魅せることができるぶんインパクトがあります.

用途は限られますが,プレゼンテーション等の演出に使えるかもしれません.

HTML TABLEの可視化

ウェブ上にはTableタグを用いて公開されるデータも多く, Tableタグを処理できるようになると,より幅広くオープンデータが利用できるようになります.

RのXMLパッケージにはそのものずばりの readHTMLTable 関数があります. これを用いると,HTMLのTABLEタグをdata.frameとして読み込むことができます.

それでは具体例をもとにreadHTMLTable関数を使ってみましょう.

東京都の平均気温の可視化

気象庁のオープンデータを利用して東京都の平均気温の傾向を可視化します.

「過去の気象データの検索」にアクセスします.

http://www.data.jma.go.jp/obd/stats/etrn/index.php

「都道府県・地方を選択」をクリックします.

「東京」をクリックします.

「地点の選択」でも「東京」をクリックします.

「観測開始から月ごとの値を表示」をクリックします.

「一覧表」「日平均気温」が選択されていることを確認します.

このTableタグをdata.frameとして読み込み、可視化します.

if(!require(reshape2)){
  install.packages("reshape2")
  library(reshape2)
}
if(!require(XML)){
  install.packages("XML")
  library(XML)
}
if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
# HTMLファイルをダウンロードします
if(!file.exists("data00.html")){
  download.file(url="http://www.data.jma.go.jp/obd/stats/etrn/view/monthly_s3.php?prec_no=44&block_no=47662&year=2015&month=&day=&view=a1",destfile="data00.html",method = "auto")
}
# HTMLページからTableデータを取り出す
tmp <- XML::readHTMLTable(readLines("data00.html"))
# 処理対象のTableを取り出す(HTMLページに複数のTableがある場合はtablefix1,tablefix2...のようになる)
tbl <- tmp$tablefix1
#列名を設定する
colnames(tbl)<-c("Year","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Ave")
data <- reshape2::melt(tbl, id=c("Year", "Ave"))
# Factor型を数値に変換する
data$Ave<-as.numeric(as.character(data$Ave))
data$value<-as.numeric(as.character(data$value))
data$Year<-as.numeric(as.character(data$Year))

g<-ggplot2::ggplot(data=data, ggplot2::aes(x = Year, y = value,group=variable)) + geom_path()
# 月ごとにグラフを表示する
g <- g + ggplot2::facet_wrap(~variable )
# X軸目盛りの範囲を1875-2015までとし、50年ごとに目盛りを入れる
g <- g +ggplot2::scale_x_continuous(breaks=seq(1875,2015,by=50),limits=c(1875,2015))
# 回帰直線を描画する
g <- g + ggplot2::stat_smooth(method = "lm")
# 実際の描画
plot(g)

以上のように readHTMLTable 関数を用いると,Tableタグから簡単にデータを読み取ることができます.

ところで,可視化結果を見ると,東京都の平気気温は徐々に上昇しているようです.短期的(といっても100年程度)には温暖化傾向にあると考えることもできそうです.

おわりに

本日の演習では Google Earth を使った可視化,HTML Table からのデータ読み込み&可視化を体験しました. 駆け足でしたがRでのオープンデータ処理を体験してもらいました. ウェブ上に存在する多くのオープンデータを処理可能な手法をお伝えしました. 消化できていないことも多いと思いますが,ヘルプを参照する(?geom_path or ??geom_path などして参照可能 ),サンプルプログラムを改造する, ウェブリソースを参照する,などして是非とも理解を深めてみてください.

次回からRをすこし離れ,ウェブアプリを使ったオープンデータの可視化手法に挑戦します. 今回の課題は自由度が高いため,難しいかと思いますが,がんばってみてください. おもしろい可視化を楽しみにしています.

課題

本日の課題はRを用いたオープンデータの可視化です.これまで取り上げたオープンデータ以外であれば何でも構いません.これまで取り上げたオープンデータであっても,可視化手法が異なる場合はOKとします(できれば新しいデータで挑戦してみてください).

提出物は,可視化に用いたRスクリプト,利用したオープンデータ・可視化結果についての説明文です.

以降はアイデアのヒントにしてみてくでさい.

様々なデータの可視化

風向きの可視化(東京都)

オープンデータの可視化例として「東京風速(http://air.nullschool.net/)」があります.

こういったお洒落な可視化は手間がかかりますが,頑張れば近いものは実現できます.

今回は「東京風速」のアイデアをもとに,東京都における風速と風向きの可視化を行います.

東京都Shapeファイルの描画

esriジャパンから「全国市区町村界データ」(shapefile)をダウンロードします.

http://www.esrij.com/products/japan-shp/

このShapefileには市区町村コード(JCODE)が含まれています.このコードを利用して東京都のShapefileを取り出します.

JCODEについては総務省:全国地方公共団体コード(http://www.soumu.go.jp/denshijiti/code.html)で確認することができます.

ダウンロード後,Rから読み込み可能な場所にShapefileを展開します.

if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
if(!require(XML)){
  install.packages("XML")
  library(XML)
}
if(!require(jsonlite)){
  install.packages("jsonlite")
  library(jsonlite)
}
if(!require(broom)){
  install.packages("broom")
  library(broom)
}
if(!require(dplyr)){
  install.packages("dplyr")
  library(dplyr)
}
shp<-maptools::readShapePoly("japan_ver80.shp")
#市区町村コード(JCODE)を数値に変換します
shp@data$JCODE<-as.numeric(as.character(shp@data$JCODE))
# 文字コードを変換します
shp@data$KEN<-iconv(shp@data$KEN,from="SJIS")
shp@data$SIKUCHOSON<-iconv(shp@data$SIKUCHOSON,from="SJIS")
shp@data$SICHO<-iconv(shp@data$SICHO,from="SJIS")
shp@data$GUN<-iconv(shp@data$GUN,from="SJIS")
shp@data$SEIREI<-iconv(shp@data$SEIREI,from="SJIS")
# 東京都(三宅島等を除く)の領域データだけを取り出します.
tokyo<-subset(shp,shp@data$JCODE>13100&shp@data$JCODE<13361)
points <- broom::tidy(tokyo)
tokyo@data$id <- rownames(tokyo@data)
tokyo.df <- dplyr::inner_join(points, tokyo@data,by = "id")
ggplot2::ggplot(tokyo.df)+ggplot2::geom_polygon(ggplot2::aes(long,lat,group=group))

観測地点の可視化

観測地点の位置情報を取得します.

ちょっとずるいですが,ここでは東京風速のソースコード(https://github.com/cambecc/air)から観測地点の位置情報を取得します.

#「東京風速」から観測地点が記録されたjsonファイルを取得する
if(!file.exists("station-data.json")){
  download.file(url="https://raw.githubusercontent.com/cambecc/air/master/station-data.json",destfile="station-data.json",method = "auto")
}
# JSONファイルを読み込む
station<-as.data.frame(jsonlite::fromJSON("station-data.json"))
colnames(station)<-c("id","name","addr","lat","lon")
# Factorを数値に反韓する
station$lon <- as.numeric(as.character(station$lon))
station$lat <- as.numeric(as.character(station$lat))
station$id <- as.numeric(as.character(station$id))
# 観測地点を描画する
ggplot2::ggplot()+ggplot2::geom_polygon(data=tokyo.df,ggplot2::aes(long,lat,group=group))+ggplot2::geom_point(data=station,ggplot2::aes(lon,lat),colour="white")

風速と風向きの可視化

「大気汚染地図情報(速報値)(http://www.taiki.kankyo.metro.tokyo.jp/)」から風速と風向きを取得します.

このページのデータは基本Tableタグで提供されていますので,本日のテクニック(readHTMLTable関数)が利用できます.

#UNIXTIMEを取得する
unixtime <- trunc(as.numeric(Sys.time()))
#一般局のhtmlページを取得する
#URLを作成する
url <-paste("http://www.taiki.kankyo.metro.tokyo.jp/cgi-bin/bunpu1/p160.cgi?wd=",unixtime,"=",unixtime,"=1=800=2====2=",sep="")
if(!file.exists("taiki_tokyo1.html")){
  download.file(url=url,destfile="taiki_tokyo1.html",method = "auto")
}
tbl1<-XML::readHTMLTable("taiki_tokyo1.html",header = FALSE)
#自排局のhtmlページを取得する
url <-paste("http://www.taiki.kankyo.metro.tokyo.jp/cgi-bin/bunpu1/p160.cgi?wd=",unixtime,"=",unixtime,"=2=800=2====2=",sep="")
if(!file.exists("taiki_tokyo2.html")){
  download.file(url=url,destfile="taiki_tokyo2.html",method = "auto")
}
tbl2<-XML::readHTMLTable("taiki_tokyo2.html",header = FALSE)
#縦方向にTableを連結する
tbl<-rbind(data.frame(tbl1[4]),data.frame(tbl2[4]))
#列名を設定する
colnames(tbl)<-c("no","id","name","SO2","OX","NO","NO2","NOx","CO","SPM","NMHC","CH4","PM25","WD","WS","TMP","HUM","INS")
tbl$id <- as.numeric(as.character(tbl$id))
#観測地点のデータとidをキーに結合する
tmp<-dplyr::inner_join(station,tbl,by="id")

#風向きを角度に変換するためのdata.frameを用意する
lbl<-c("&nbsp","C","E","ENE","NE","NNE","N","NNW","NW","WNW","W","WSW","SW","SSW","S","SSE","SE","ESE")
#対応する角度
deg<-c(-1,-1,seq(0,337.5,22.5))
wd<-data.frame(lbl,deg)

#データの風向きを角度に置き換える
levels(tmp$WD)<-wd$deg
tmp$WD<-as.numeric(as.character(tmp$WD))
tmp$WS<-as.numeric(as.character(tmp$WS))
#Naは0にする
tmp$WS[is.na(tmp$WS)]<-0

#風向きと風速から矢印を描画するためのデータ生成関数
arrow_data<-function(sx,sy,len,deg){
  # x = len * cos(deg*pi/180), y = len * sin(deg*pi/180)
  return(data.frame(x=c(sx,len*cos(deg*pi/180)+sx),y=c(sy,len*sin(deg*pi/180)+sy)))
}
#矢印を描画する(x=0,y=0,風速10,角度45)
ggplot2::ggplot(arrow_data(0,0,10,45))+ggplot2::geom_path(aes(x,y),arrow=grid::arrow())

#観測地点を起点に風速と風向きから矢印描画に必要なデータを生成する
tmp<-cbind(tmp,arrow_data(tmp$lat,tmp$lon,tmp$WS/1000,tmp$WD))
#風が吹いたデータのみ抽出する
w<-subset(tmp,tmp$WD!=-1)
#描画に必要な列だけ取り出す
w2<-w[c(1,4,5,23,24)]
#縦方向にデータを連結する
w3<-rbind(as.matrix(w2[c(1,2,3)]),as.matrix(w2[c(1,4,5)]))
rownames(w3)<-seq(1:nrow(w3))
w3<-as.data.frame(w3)
colnames(w3)<-c("id","lat","lon")

g<-ggplot2::ggplot()+ggplot2::geom_polygon(data=tokyo.df,ggplot2::aes(long,lat,group=group))
# 矢印は緑色
g<-g+ggplot2::geom_path(data=w3,ggplot2::aes(x=lon,y=lat,group=id),arrow=grid::arrow(length = grid::unit(5,units = "point")),colour="green")
g<-g+ggplot2::theme_bw()
plot(g)

「東京風速」のようにはいきませんが,風向きの可視化ができました. htmlファイルを直接読み込む形にすれば,実行するたびにその時点の風向きを可視化するようになります.

台風進路の可視化

気象庁の台風進路のベストトラックデータを可視化してみます.

http://www.jma.go.jp/jma/jma-eng/jma-center/rsmc-hp-pub-eg/besttrack.html

ベストトラックのデータフォーマットは以下のとおりです.

http://www.jma.go.jp/jma/jma-eng/jma-center/rsmc-hp-pub-eg/Besttracks/e_format_bst.html

フォーマットの日本語訳は以下のページで参照できます.

http://homepage3.nifty.com/typhoon21/general/bst-format.html

ベストトラックはヘッダー行とデータ行,二種類のフォーマットが混在するデータ形式です.

まず,ベストトラックデータをダウンロードします.

if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
if(!require(maps)){
  install.packages("maps")
  library(maps)
}
if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
if(!require(mapdata)){
  install.packages("mapdata")
  library(mapdata)
}
if(!require(scales)){
  install.packages("scales")
  library(scales)
}
if(!file.exists("bst_all.zip")){
  download.file(url="http://www.jma.go.jp/jma/jma-eng/jma-center/rsmc-hp-pub-eg/Besttracks/bst_all.zip",destfile="bst_all.zip",method = "auto")
  unzip("bst_all.zip",exdir = ".")
} 

ヘッダー行を取り出します.read.fwf関数で固定長データとして読み込みます.

#一行ずつ読み込む
tmp <- readLines("bst_all.txt")
#66666で始まるヘッダー行のみを取り出す
header_line_no <- grep("^66666", tmp)
#Header: 1:Indicator, 2:International number ID, 3:Number of data lines, 4:Tropical cyclone number ID, 5:International number ID, 6:Flag of the last data line, 7:Difference between the time of the last data and the time of the final analysis, 8:Name of the storm, 9:Date of the latest revision
header<-read.fwf(textConnection(tmp[header_line_no]),widths = c(6,5,4,5,5,2,2,21,9),col.names = c("IND","ID","N","TCID","ID2","FLG","DEF","NAME","REV"),fill=TRUE,header = FALSE,colClasses="character")

データ行を取り出します.

#entry: 1:Time of analysis, 2:Grade, V3:Latitude of the center, 4:Longitude of the center, 5:Central pressure, 6:Maximum sustained wind speed, 7:Direction of the longest radius of 50kt winds or greater, 8:The longest radius of 50kt winds or greater, 9:The shortest radius of 50kt winds or greater, 10:Direction of the longest radius of 30kt winds or greater, 11:The longest radius of 30kt winds or greater, 12:The shortest radius of 30kt winds or greater, 13:Indicator of landfall or passage
entry <- read.fwf(textConnection(tmp[-header_line_no]),widths=c(9,4,2,4,5,5,4,2,5,5,2,5,5,2),col.names = c("TIME","IND","GRADE","LAT","LON","hPa","MWS","WSD","WSLR","WSSR","SWD","SWLR","SWSR","SIZE","STR","LND"),header = FALSE,fill=TRUE,colClasses="character")

ヘッダーとデータをつなぎ合わせます.

#ヘッダーにIDをつけます
entry$ID<-rep(header$ID,header$N)
#IDをキーにヘッダーとデータをマージします
data<-merge(header,entry,by= "ID")
#観測年(西暦2桁)を西暦4桁に変換します
#先頭二文字を取り出します
data$YEAR<-as.numeric(substring(data$TIME,0,2))
#15以上の場合は1900を足します
data$YEAR[data$YEAR>16]<-data$YEAR[data$YEAR>15]+1900
#15以下の場合は2000を足します
data$YEAR[data$YEAR<=16]<-data$YEAR[data$YEAR<=15]+2000
#日時データとして処理できるようにします
data$TIME[data$YEAR>=2000]<-paste("20",data$TIME[data$YEAR>=2000],sep = "")
data$TIME[data$YEAR<2000]<-paste("19",data$TIME[data$YEAR<2000],sep = "")
data$TIME<-strptime(data$TIME, "%Y%m%d%H")
#緯度、軽度を修正します
data$LAT <- as.numeric(data$LAT)/10
data$LON <- as.numeric(data$LON)/10
#台風名に両端の空白を削除します
# trimws はR 3.2.0以降で利用可能
#data$NAME <- trimws(data$NAME,which = "both")
data$hPa <- as.numeric(data$hPa)

#日本周辺地図を取り出します
LON_LIMIT <- scales::expand_range(range(data$LON), add=0)
LAT_LIMIT <- scales::expand_range(range(data$LAT), add=0)
jpn <- data.frame(maps::map(plot=FALSE, xlim = LON_LIMIT, ylim = LAT_LIMIT)[c("x","y")])
ggplot2::ggplot()+ggplot2::geom_path(data=jpn,ggplot2::aes(x, y, colour = NULL)) 

2016年

#2016-01-01から2016-12-31までのデータを取り出す
data2<-subset(data,data$TIME>as.POSIXlt("2016-01-01")&data$TIME<as.POSIXlt("2016-12-31"))
p <- ggplot2::ggplot(data2, ggplot2::aes(LON, LAT, colour = NAME)) + 
  ggplot2::geom_point(ggplot2::aes(size = hPa,fill=NAME), shape = 21, alpha = 0.2) +
  ggplot2::geom_path() +
  ggplot2::geom_path(ggplot2::aes(x, y, colour = NULL), jpn) +
  ggplot2::coord_map() +
  ggplot2::scale_size_continuous(range= c(1, 10), trans = "reverse") +
  ggplot2::theme_bw()+ggplot2::theme(legend.position="none")+ggplot2::ggtitle("2016")
plot(p)

2015年

#2015-01-01から2015-12-31までのデータを取り出す
data2<-subset(data,data$TIME>as.POSIXlt("2015-01-01")&data$TIME<as.POSIXlt("2015-12-31"))
p <- ggplot2::ggplot(data2, ggplot2::aes(LON, LAT, colour = NAME)) + 
  ggplot2::geom_point(ggplot2::aes(size = hPa,fill=NAME), shape = 21, alpha = 0.2) +
  ggplot2::geom_path() +
  ggplot2::geom_path(ggplot2::aes(x, y, colour = NULL), jpn) +
  ggplot2::coord_map() +
  ggplot2::scale_size_continuous(range= c(1, 10), trans = "reverse") +
  ggplot2::theme_bw()+ggplot2::theme(legend.position="none")+ggplot2::ggtitle("2015")
plot(p)

日本地図だけで描画する場合(要mapdataライブラリ)

jpn <- data.frame(maps::map("japan",plot = FALSE)[c("x","y")])
p <- ggplot2::ggplot(data2, ggplot2::aes(LON, LAT, colour = NAME)) + 
  ggplot2::geom_point(ggplot2::aes(size = hPa,fill=NAME), shape = 21, alpha = 0.2) +
  ggplot2::geom_path() +
  ggplot2::geom_path(ggplot2::aes(x, y, colour = NULL), jpn) +
  ggplot2::coord_map() +
  ggplot2::scale_size_continuous(range= c(1, 10), trans = "reverse") +
  ggplot2::theme_bw()+ggplot2::theme(legend.position="none")
plot(p)

2015年は関東を直撃した台風はほとんどなかったようです.

インフルエンザ発症数の可視化(川崎市)

今年のインフルエンザはA型、B型どちらが流行るのでしょうか.

オープンデータを用いて傾向をみてみましょう.

川崎市感染症情報発信システムにアクセスします.

https://kidss.city.kawasaki.jp/modules/topics/

[リアルタイムサーベイランス]>[ダウンロード]の順番をクリックします.

ダウンロードサービスの各項目を以下のように設定します.

  • 平成26年(2014年)1月1日から平成28年(2016年)12月31日まで
  • 集計方法:「年齢階級別集計」
  • データの種類:「日別データ」

疾患選択で「A型インフルエンザ」を選び、CSVダウンロードをクリックします.

疾患選択で「B型インフルエンザ」を選び、CSVダウンロードをクリックします.

ダウンロードしたCSVファイルをRから参照可能なディレクトリに保存します. ファイル名を変更しておくとあとあと便利です.(下記ではkawasaki_fluA.csv,kawasaki_fluB.csvとしている)

if(!require(dplyr)){
  install.packages("dplyr")
  library(dplyr)
}
if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
#文字コードはSJIS、最初の2行はスキップする
A<-read.csv("kawasaki_fluA.csv",fileEncoding = "SJIS",header = TRUE,skip = 2)
B<-read.csv("kawasaki_fluB.csv",fileEncoding = "SJIS",header = TRUE,skip = 2)
colnames(A)<-c("year","month","day","m5","m11","y1","y2","y3","y4","y5","y6","y7","y8","y9","y14","y19","y29","y39","y49","y59","y69","y79","y80over","total")
colnames(B)<-c("year","month","day","m5","m11","y1","y2","y3","y4","y5","y6","y7","y8","y9","y14","y19","y29","y39","y49","y59","y69","y79","y80over","total")
A$type<-"A"
B$type<-"B"
d<-rbind(A,B)
d$total <- as.numeric(as.character(d$total))

#データをグルーピング、合計を計算して、可視化
dat <- d%>%dplyr::group_by(year,month,type)%>%dplyr::summarise(total=sum(total))
ggplot2::ggplot(data = dat,mapping = ggplot2::aes(x=month,y=total,group=type,colour=type))+ggplot2::geom_path()+ggplot2::facet_wrap(~year)+ggplot2::scale_x_continuous(breaks=seq(1,12,by=1))+ggplot2::ggtitle("川崎市のインフルエンザ患者数")+ggplot2::theme_gray(base_family = "yugo")

今年はA型?まだ爆発的とはいっていないようです.

青空文庫:夏目漱石「こころ」の可視化

青空文庫というウェブサイトがあります.

http://www.aozora.gr.jp/

青空文庫は、誰にでもアクセスできる自由な電子本を、図書館のようにインターネット上に集めようとする活動です。著作権の消滅した作品と、「自由に読んでもらってかまわない」とされたものを、テキストとXHTML(一部はHTML)形式に電子化した上で揃えています。

ref: http://www.aozora.gr.jp/guide/aozora_bunko_hayawakari.html

今回は青空文庫から夏目漱石の「こころ」のテキストデータをダウンロードし,可視化を試みます.

http://www.aozora.gr.jp/cards/000148/files/773_ruby_5968.zip

なお,テキストデータは形態素解析ソフトウェア Mecab(http://taku910.github.io/mecab/)で処理します. 本演習では,Mecabはインストール済みであることを想定しています.

「773_ruby_5968.zip」を展開し,「kokoro.txt」を形態素解析します.

コマンドプロンプトで以下を実行してください.

mecab kokoro.txt > output.txt

「output.txt」には形態素解析結果が出力されます. この処理結果をRで可視化します.

if(!require(wordcloud)){
  install.packages("wordcloud")
  library(wordcloud)
}
tbl<-read.table("output.txt",header = FALSE,fileEncoding = "SJIS",fill = TRUE)
# 品詞:名詞で一般名詞とラベル付けされている語を取り出す
verb <- tbl[grep("^名詞,一般",tbl$V2),]
# tableで語の出現頻度を計算し,data.frameに変換します.
df <- as.data.frame(table(verb$V1))
# データを降順に並び替え,上位30の語を取り出します
df2<-df[order(df$Freq,decreasing = TRUE),][1:30,]
# ワードクラウド表示します
wordcloud::wordcloud(df2$Var1,df2$Freq,scale = c(6, .2),random.order = T, rot.per = .15, colors = rainbow(5))

なんとなく作品のテーマがあらわれていますね.ちなみに,こうした可視化手法をワードクラウドといいます.

パナマ文書:ネットワークデータの可視化

パナマ文書についての可視化にも挑戦してみます(厳密にはパナマ文書を含む).

※データサイズが大きいので注意してください.

ICIJ Offshore Leaks Database(https://offshoreleaks.icij.org/)にアクセスします.

メニュにある「DOWNLOAD」をクリックします.

「How to download this database」の第二パラグラフ最後の「 You may download an archive of all these files here.」の「here」をクリックし, データベース「data-csv.zip」をダウンロードします.

ダウンロードしたファイルを展開します.

データベースは以下の5つのファイルから構成されます.

  • Addresses.csv : 住所リスト
  • Entities.csv : 法人リスト
  • Intermediaries.csv : 仲介者リスト
  • Officers.csv : 役員リスト
  • all_edges.csv : 上記の関係情報

今回は特定の人物にフォーカスして ネットワーク図を描画してみたいと思います. ネットワーク図の描画にはネットワーク分析ライブラリ「igraph(http://igraph.org/)」を使います.

以下にコード例を示します.ファイルパスは適宜読み替えて下さい.

if(!require(igraph)){
  install.packages("igraph")
  library(igraph)
}
if(!require(dplyr)){
  install.packages("dplyr")
  library(dplyr)
}
#役員リストを読み込む
officers <- read.csv("data-csv/Officers.csv",quote = "\"",sep = ",",stringsAsFactors = FALSE)
#法人リストを読み込む
entities = read.csv("data-csv/Entities.csv",quote = "\"",sep = ",",stringsAsFactors = FALSE)
#住所リストを読み込む
addresses = read.csv("data-csv/Addresses.csv",quote = "\"",sep = ",",stringsAsFactors = FALSE)
#仲介者リストを読み込む
intermediaries = read.csv("data-csv/Intermediaries.csv",quote = "\"",sep = ",",stringsAsFactors = FALSE)
#上記リストの関係情報
all_edges = read.csv("data-csv/all_edges.csv",quote = "\"",sep = ",",stringsAsFactors = FALSE)

# 役員名「Sigmundur David Gunnlaugsson」のデータを取り出す
target <- officers[officers$name =="Sigmundur David Gunnlaugsson",]

# 上記ターゲットと関連あるnode_idを取り出す
edges <- all_edges %>% dplyr::filter(node_1 %in% target$node_id | node_2 %in% target$node_id)

#探索する最大の深さ
depth <- 3
n <- 0
for(i in 1:depth){
  # 関連あるnode_idに紐づくnode_idを探す
  tmp_edges <- all_edges %>% dplyr::filter(node_1 %in% edges$node_1| node_2 %in% edges$node_1|node_2 %in% edges$node_1|node_2 %in% edges$node_2)
  # 行数を取得する
  tmp <- nrow(tmp_edges)
  edges <- rbind(edges,tmp_edges)
  # 行数に変更がない(=つながりを持つデータが存在しない)場合は繰返し処理を修了する
  if( n == tmp) break
  # 現在の行数を更新する
  n <- tmp
}


# 重複データを削除する
edges <- unique(edges)
nrow(edges)
## [1] 25
# ノードを抽出する
vertices<-c(edges$node_1,edges$node_2)
# 重複データを削除する
vertices<-unique(vertices)
# data.frameに変換する
vertices<-data.frame(node_id=vertices)

#node_idをキーにデータを統合する(name列だけ).
tmp<-dplyr::inner_join(vertices,officers[,c('node_id','name')],by='node_id')
tmp$color <- "chartreuse"
tmp$shape<-"circle"

tmp2<-dplyr::inner_join(vertices,entities[,c('node_id','name')],by='node_id')
tmp2$color <- "Pink"
#tmp2$shape<-"rectangle"
tmp2$shape<-"circle"

tmp3<-dplyr::inner_join(vertices,intermediaries[,c('node_id','name')],by='node_id')
tmp3$color <- "Lightgreen"
#tmp3$shape <- "csquare"
tmp3$shape <- "circle"

tmp4<-dplyr::inner_join(vertices,addresses[,c('node_id','address')],by='node_id')
colnames(tmp4)<-c('node_id','name')
tmp4$color <- "cadetblue"
#tmp4$shape <- "vrectangle"
tmp4$shape <- "circle"


#データフレームを統合する
vertices2<-rbind(tmp,tmp2,tmp3,tmp4)

#ネットワークグラフ作成のためにエッジデータを最定義する
edges2<-data.frame(from=edges$node_1,to=edges$node_2,rel=edges$rel_type)

# 無向グラフデータを作成する
g<-igraph::graph.data.frame(d=edges2,directed = FALSE,vertices = vertices2)

#各ノードの次数(他のノードとどれだけ繋がっているか)を取得する
g.ret <- igraph::degree(g)

#グラフを描画する
igraph::plot.igraph(g,
     # 次数の5倍をノードサイズとする
     vertex.size=g.ret*5,
     # ノードの形はshape要素で指定する
     vertex.shape=V(g)$shape,
     # ノードのラベルはname要素で指定する
     vertex.label=V(g)$name,
     # ノードの色はcolor要素で指定する
     vertex.color=V(g)$color,
     # ラベルの色は「gray50」とする
     vertex.label.color="gray50",
     # ノードの枠は「white」とする
     vertex.frame.color="white",
     # ノードラベルのサイズを設定する
     vertex.label.cex=0.5,
     # エッジの色は「gray80」とする
     edge.color="gray80")

グラフ描画の際「tkplot」を使うと, インタラクティブにグラフを操作することができます.

igraph::tkplot(g,
     # 次数の5倍をノードサイズとする
     vertex.size=g.ret*5,
     # ノードの形はshape要素で指定する
     vertex.shape=V(g)$shape,
     # ノードのラベルはname要素で指定する
     vertex.label=V(g)$name,
     # ノードの色はcolor要素で指定する
     vertex.color=V(g)$color,
     # ラベルの色は「gray50」とする
     vertex.label.color="gray50",
     # ノードの枠は「white」とする
     vertex.frame.color="white",
     # ノードラベルのサイズを設定する
     vertex.label.cex=0.5,
     # エッジの色は「gray80」とする
     edge.color="gray80")

パナマ文書により辞任に追い込まれた 元アイスランド首相「Sigmundur David Gunnlaugsson」に関するデータです. これにより何が言えるというわけでもありませんが, ネットワーク解析を使うと,人と人など関係性をもったデータの分析ができるようになります.

興味がある方はigraphのチュートリアル的なデモをみてみてください.

igraph::igraphdemo("community")

番外編:Linked Open Data(LOD)

演習では 5 Star OPEN DATA の4段目と5段目については取り上げませんでしたが, 4段目以降はより高度な処理を目指した意欲的な内容となっています.

簡単に表現すると,4段目のRDF(Resource Description Framework)でデータを表現することを提案しています.RDFはウェブにおけるデータの表現モデルの一つで, データをURI(Uniform Resource Identifier)で識別し,そのURIとURIの関係性によりデータを表現します.

5段目はデータをLOD(Linked Open Data)として表現することを提案しています.LODはRDFデータとRDFデータの関係性を記述し,データのネットワークを構築しようとする意欲的な試みとなります.

興味がある方はググったり,下記サイトを参照してみて下さい.

ここでは Linked Open Data な世界を体験するために, RでSPARQL(SPARQL Protocol and RDF Query Language)を実行してみます.

まだまだ、LODによるオープンデータ公開事例は少ないですが、その可能性を感じてみてください.

dbpedia.org

wikipediaのLOD化を推進する dbpedia.org(http://wiki.dbpedia.org/) に対するSPARQLクエリです。

日本生まれ

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://dbpedia.org/sparql"
q<-"SELECT *
WHERE {
?s <http://dbpedia.org/ontology/birthPlace> <http://dbpedia.org/resource/Japan> .
}"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

1980年以降に誕生した日本生まれの人

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://dbpedia.org/sparql"
q<-"SELECT *
WHERE {
?s <http://dbpedia.org/ontology/birthPlace> <http://dbpedia.org/resource/Japan> .
?s <http://dbpedia.org/ontology/birthDate> ?date .
filter (?date >= '1980-1-1'^^<http://www.w3.org/2001/XMLSchema#date>) .
}"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

3000m以上の山

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://dbpedia.org/sparql"
q<-"SELECT * WHERE {
?s <http://dbpedia.org/ontology/elevation> ?elevation .
?s <http://dbpedia.org/property/wordnet_type> <http://www.w3.org/2006/03/wn/wn20/instances/synset-mountain-noun-1> .
?s <http://dbpedia.org/property/name> ?name .
filter (?elevation >= 3000 )
}
order by desc(?elevation)"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

WikipediaのページURLからDBPediaのリソースURI

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://dbpedia.org/sparql"
q<-"SELECT * WHERE {
?s <http://xmlns.com/foaf/0.1/isPrimaryTopicOf> <http://en.wikipedia.org/wiki/Japan>
}"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

日本の都道府県、人口、緯度経度

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://dbpedia.org/sparql"
q<-"SELECT DISTINCT ?title ?lat ?lon ?populationTotal where{
?s <http://www.w3.org/1999/02/22-rdf-syntax-ns#type> <http://dbpedia.org/class/yago/PrefecturesOfJapan> .
?s <http://www.w3.org/2000/01/rdf-schema#label> ?title .
?s <http://dbpedia.org/ontology/populationTotal> ?populationTotal .
?s <http://www.w3.org/2003/01/geo/wgs84_pos#lat> ?lat .
?s <http://www.w3.org/2003/01/geo/wgs84_pos#long> ?lon .
filter langMatches(lang(?title),'en') .
}"

res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

ja.dbpedia.org

dbpedia日本語版(http://ja.dbpedia.org/)では日本語を使ったSPARQL検索ができます。

東京都

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://ja.dbpedia.org/sparql" 
q <- "SELECT DISTINCT * WHERE 
{ 
<http://ja.dbpedia.org/resource/東京都> ?p ?o .
 }"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

ロック音楽のリスト(もしあれば画像uriも)

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://ja.dbpedia.org/sparql" 
q <- "SELECT DISTINCT ?label ?depiction
WHERE {
  ?s <http://ja.dbpedia.org/property/genre> <http://ja.dbpedia.org/resource/ロック_(音楽)> .
  ?s <http://www.w3.org/2000/01/rdf-schema#label> ?label .
  OPTIONAL { ?s <http://xmlns.com/foaf/0.1/depiction> ?depiction .  }
}"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

誕生日が1月1日

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://ja.dbpedia.org/sparql" 
q <- "SELECT DISTINCT ?label ?birthYear
WHERE {
?s <http://ja.dbpedia.org/property/生日> 1 .
?s  <http://ja.dbpedia.org/property/生月> 1 .
?s   <http://ja.dbpedia.org/property/生年> ?birthYear .
?s    <http://www.w3.org/2000/01/rdf-schema#label> ?label .
}
ORDER BY ?birthYear"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

全国の地域限定ゆるキャラ

if(!require(SPARQL)){
  install.packages("SPARQL")
  library(SPARQL)
}
endpoint <- "http://ja.dbpedia.org/sparql" 
q <- "SELECT DISTINCT ?mascot ?name
WHERE {
?areamascots <http://www.w3.org/2004/02/skos/core#broader> <http://ja.dbpedia.org/resource/Category:地域限定のマスコット> .
?areamascots  <http://www.w3.org/2000/01/rdf-schema#label> ?areaname .
?mascot <http://dbpedia.org/ontology/wikiPageWikiLink> ?areamascots .
?mascot  <http://www.w3.org/2000/01/rdf-schema#label> ?name.
}"
res <- SPARQL(url=endpoint,query=iconv(q,to="utf8"))
head(res)

戻る


クリエイティブ・コモンズ・ライセンス
Masaharu Hayashi を著作者とするこの 作品 は クリエイティブ・コモンズの 表示 4.0 国際 ライセンスで提供されています。